#include "HighOrderApproximation.h"

namespace PRS{

	HighOrderApproximation::HighOrderApproximation(){
	}

	HighOrderApproximation::HighOrderApproximation(SimulatorParameters* pSimPar){
		pNodeSL = getNodeSlopeLimitFunc(pSimPar);
		pEdgeSL = getEdgeSlopeLimitFunc(pSimPar);
		koef = pSimPar->get_koef();
	}

	HighOrderApproximation::~HighOrderApproximation(){
	}

	// Initialize a pointer to node slope limiter
	NodeSlopeLimiter* HighOrderApproximation::getNodeSlopeLimitFunc(SimulatorParameters* pSimPar){
		int val = pSimPar->getNodeSlopeLimitFunc();
		if (val==node_MUSCL || val==node_Superbee || val== node_Minmod || val==node_Osher || val==node_Van_Albada){
			return new Node_MUSCL();
		}
//		else if (val==node_WoodField){
//			return new Node_WoodField();
//		}
		else{
			throw Exception(__LINE__,__FILE__,"Unknown high order methods.\n");
		}
	}

	EdgeSlopeLimiter* HighOrderApproximation::getEdgeSlopeLimitFunc(SimulatorParameters* pSimPar){
		int val = pSimPar->getNodeSlopeLimitFunc();
		if (val==MUSCL || val==SUPERBEE || val== MINMOD || val==OSHER || val==VAN_ALBADA){
			return new Edge_MUSCL();
		}
//		else if (val==WOODFIELD){
//			return new Edge_WoodField();
//		}
		else{
			throw Exception(__LINE__,__FILE__,"Unknown high order methods.\n");
		}
	}

	double HighOrderApproximation::getSw_HighOrderApproximation(SimulatorParameters* pSimPar, EdgeData* edata, int dim, double &Sw_I, double &Sw_J){
		static ofstream fid;
//		static bool key=true;
//		if(key){
//			fid.open("HO.txt");
//			key=false;
//		}


		// start CPU time clock
		double startt = MPI_Wtime();

		double delta_Sw_JI = Sw_J - Sw_I;

		double SLII,SLJJ;			// edge slope limiter
		double slimit_I, slimit_J;	// nodal slope limiter
		double DSwII, DSwJJ;

		pNodeSL->calculateNodeSlopeLimiters(edata->vertex_0,edata->vertex_1,slimit_I,slimit_J);
		pEdgeSL->calculateEdgeSlopeLimiters(pSimPar,edata,dim,Sw_I,Sw_J,SLII,SLJJ,DSwII,DSwJJ);

		//fid << slimit_I << "\t" << slimit_J << "\t" << Sw_I << "\t" << Sw_J <<"\t"<< SLII << "\t" << SLJJ << Sw_J <<"\t"<< DSwII << "\t" << DSwJJ << endl;
		//fid << edata->vertex_0->ID << "\t" << edata->vertex_1->ID << "\t" << SLII << "\t" << SLJJ <<"\t"<< DSwII << "\t" << DSwJJ << endl;

		// final limiter function
		const double SLI = SLII*slimit_I;
		const double SLJ = SLJJ*slimit_J;

		/*
		 * Modified Taylor expansion series extrapolating nodal saturation value
		 * on volume control interfaces. Injection wells are excluded.
		 */
		if ( !pSimPar->isInjectionWell( edata->vertex_0->flag ) ){
			Sw_I = Sw_I + (SLI/4.)*((1.-koef)*DSwII + (1.+koef)*delta_Sw_JI);
		}
		if ( !pSimPar->isInjectionWell( edata->vertex_1->flag ) ){
			Sw_J = Sw_J - (SLJ/4.)*((1.-koef)*DSwJJ + (1.+koef)*delta_Sw_JI);
		}
		//printf("Sw_I: %.5E->%.5E\tSw_J: %.5E->%.5E\tDSwII: %.5E\tDSwJJ: %.5E\tdelta_Sw_JI: %.5E\tSLI:%.5E\tSLJ:%.5E\n",swi,Sw_I,swj,Sw_J,DSwII,DSwJJ,delta_Sw_JI,SLI,SLJ);

		// finish CPU time clock. Return CPU time consumption (seconds)
		double endt = MPI_Wtime();
		return endt-startt;
	}
}
